##########################################################################################

library(data.table)
library(optparse)
library(dplyr)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--sample_public_info"), type = "character") ,
    make_option(c("--cancer_file"), type = "character") ,
    make_option(c("--cancer_nmu_file"), type = "character") ,
    make_option(c("--sample_info"), type = "character") ,
    make_option(c("--gtf_file"), type = "character") ,
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    gene <- "MUC6"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    cancer_file <- paste(work_dir,"/maf_public/All_use.addVAF.maf",sep="")
    cancer_nmu_file <- paste(work_dir,"/maf/All_GGA.cancer.maf",sep="")
    sample_public_info <- paste(work_dir,"/public_ref/combine/MutationInfo.combine.addMolecularSubType.tsv",sep="")
    sample_info <- "~/20220915_gastric_multiple/dna_combinePublic/baseTable/STAD_Info.addBurden.MSI_MSS.addCNVType.tsv"
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/lollipop" , "/" , gene)
    gtf_file <- "~/ref/GTF/gencode.v19.annotation.exonNum.gtf"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_public_info <- opt$sample_public_info
cancer_file <- opt$cancer_file
cancer_nmu_file <- opt$cancer_nmu_file
sample_info <- opt$sample_info
gtf_file <- opt$gtf_file
gene <- opt$gene
out_path <- opt$out_path

###########################################################################################

dir.create(out_path , recursive = T)

###########################################################################################

dat_Cancerous <- fread(cancer_file  ,sep = "\t",  quote="" ,header = T)
dat_Cancerous_nmu <- data.frame(fread(cancer_nmu_file))

dat_info <- data.frame(fread(sample_info))
dat_info_public <- data.frame(fread(sample_public_info))
dat_gtf <- fread(gtf_file)

###########################################################################################

Variant_Types <- c("Missense_Mutation","Nonsense_Mutation","Frame_Shift_Ins","Frame_Shift_Del","In_Frame_Ins","In_Frame_Del","Splice_Site","Nonstop_Mutation")

###########################################################################################
## 突变整理
dat_Cancerous$t_alt_count <- as.numeric(dat_Cancerous$t_alt_count)
dat_Cancerous$t_ref_count <- as.numeric(dat_Cancerous$t_ref_count)
dat_Cancerous <- subset(dat_Cancerous , Hugo_Symbol == gene & Variant_Classification %in% Variant_Types & ( t_alt_count > 0 | From == "TMUCIH" ) )

dat_Cancerous_nmu$t_alt_count <- as.numeric(dat_Cancerous_nmu$t_alt_count)
dat_Cancerous_nmu$t_ref_count <- as.numeric(dat_Cancerous_nmu$t_ref_count)
dat_Cancerous_nmu <- subset(dat_Cancerous_nmu , Hugo_Symbol == gene & Variant_Classification %in% Variant_Types & ( t_alt_count > 0 ) )

## 存在其它转录本
## 使用第一个转录本，因为EGFR基因的默认转录本就在第一个
maf <- data.frame(dat_Cancerous)

maf_use <- data.frame(
    gene = maf$Hugo_Symbol, 
    refseq = substring(maf$Transcript_ID , 1 , 15) ,
	chromosome = paste0( "chr" , maf$Chromosome) , start = maf$Start_position , 
    REF = maf$Reference_Allele , ALT = maf$Tumor_Seq_Allele2 ,
	aachange = maf$HGVSp_Short,
	class = maf$Variant_Classification , 
    t_alt_count = maf$t_alt_count , t_ref_count = maf$t_ref_count ,
    VAF = maf$t_alt_count/(maf$t_alt_count + maf$t_ref_count) ,
    Tumor = maf$Tumor_Sample_Barcode
    )

maf_use$refseq <- sapply(strsplit(maf_use$refseq,"[.]") , "[" , 1)

if(gene == "AK2" ){
    maf_use$refseq <- "ENST00000467905"
}

if(gene == "APC"){
    maf_use$refseq <- "ENST00000257430"
}

if( length(unique(maf_use$refseq)) > 1 & !(gene %in% c("AK2" , "APC" , "ERBB2" , "ARID1A")) ){

    maf_nmu <- dat_Cancerous_nmu
    maf_use_nmu <- data.frame(
            gene = maf_nmu$Hugo_Symbol, 
            refseq = sapply( strsplit( sapply(strsplit(maf_nmu$Other_Transcripts , "[|]") , "[" , 1) , "_") , "[" , 2 ) ,
            chromosome = paste0( "chr" , maf_nmu$Chromosome) , start = maf_nmu$Start_Position , 
            REF = maf_nmu$Tumor_Seq_Allele1 , ALT = maf_nmu$Tumor_Seq_Allele2 ,
            aachange = paste0( "p." , sapply( strsplit(sapply(strsplit(maf_nmu$Other_Transcripts , "[|]") , "[" , 1) , "p[.]") , "[" , 2 )) ,
            class = maf_nmu$Variant_Classification , 
            t_alt_count = maf_nmu$t_alt_count , t_ref_count = maf_nmu$t_ref_count ,
            VAF = maf_nmu$t_alt_count/(maf_nmu$t_alt_count + maf_nmu$t_ref_count) ,
            Tumor = maf_nmu$Tumor_Sample_Barcode
            )

    maf_use_nonnjmu <- maf_use[!maf_use$Tumor %in% maf_use_nmu$Tumor,]
    
    ## 判断所用的转录本
    refseq_use <- unique(maf_use_nonnjmu$refseq)
    index_ref <- grep( refseq_use , strsplit(maf_nmu$Other_Transcripts , "[|]")[[1]])

    maf_use_nmu <- data.frame(
            gene = maf_nmu$Hugo_Symbol, 
            refseq = sapply( strsplit( sapply(strsplit(maf_nmu$Other_Transcripts , "[|]") , "[" , index_ref) , "_") , "[" , 2 ) ,
            chromosome = paste0( "chr" , maf_nmu$Chromosome) , start = maf_nmu$Start_Position , 
            REF = maf_nmu$Tumor_Seq_Allele1 , ALT = maf_nmu$Tumor_Seq_Allele2 ,
            aachange = paste0( "p." , sapply( strsplit(sapply(strsplit(maf_nmu$Other_Transcripts , "[|]") , "[" , index_ref) , "p[.]") , "[" , 2 )) ,
            class = maf_nmu$Variant_Classification , 
            t_alt_count = maf_nmu$t_alt_count , t_ref_count = maf_nmu$t_ref_count ,
            VAF = maf_nmu$t_alt_count/(maf_nmu$t_alt_count + maf_nmu$t_ref_count) ,
            Tumor = maf_nmu$Tumor_Sample_Barcode
            )

    maf_use_nmu$refseq <- substring( maf_use_nmu$refseq , 1 , 15 )

    maf_use <- rbind( maf_use_nonnjmu , maf_use_nmu )
}

############################################################################################################
## 比较IGC和DGC的差异
dat_info_public <- subset( dat_info_public , From != "NJMU" )
dat_info_public$ID <- dat_info_public$Tumor
dat_info <- subset( dat_info , Type != "IM + IGC + DGC"& Class != "IM" )
dat_info$From <- "NJMU"
colnames(dat_info)[colnames(dat_info) == "TCGA_Class"] <- "Molecular.subtype"
colnames(dat_info)[colnames(dat_info) == "Patient"] <- "ID"

info_use <- rbind( dat_info[,c("ID" , "Tumor" , "Class" , "From" , "Molecular.subtype")] , dat_info_public[,c("ID" , "Tumor" , "Class" , "From" , "Molecular.subtype")] )

############################################################################################################
## 标记基线信息
maf_use <- merge( maf_use , info_use , by = "Tumor" )

## 若不影响蛋白质改变，则用class替换其aachange
#maf_use[maf_use$aachange=="","aachange"] <- maf_use[maf_use$aachange=="","class"]
#maf_use[maf_use$aachange=="p.NA","aachange"] <- maf_use[maf_use$aachange=="p.NA","class"]
index <- which(maf_use$class=="Splice_Site")

#maf_use[maf_use$class=="Splice_Site","aachange"] <- paste("Splice_Site" , 
#    maf_use$chromosome[index], maf_use$start[index] , maf_use$REF[index] , maf_use$ALT[index] , 
#    sep = ":")

maf_use[maf_use$class=="Splice_Site","aachange"] <- "Splice_Site" 

############################################################################################################
## 注释其外显子
## 其对应的外显子
maf_use$exon_pos <- ""
for(i in 1:nrow(maf_use)){
    exon_num <- dat_gtf[which(dat_gtf$transcript_id == maf_use$refseq[i] & dat_gtf$chr == maf_use$chromosome[i] & dat_gtf$start <= maf_use$start[i] & dat_gtf$end >= maf_use$start[i]  ),"exon_number"]
    exon_num <- as.numeric(unique(data.frame(exon_num)))
    maf_use$exon_pos[i] <- exon_num    
}

lollipop_variant <- unique(maf_use)

## 所有的信息
out_name <- paste0(out_path , "/" , gene , ".AllInfo.IGC_DGC.tsv" )
write.table(unique(lollipop_variant) , out_name , row.names=F , sep ="\t" , quote = F)

############################################################################################################
## 输出画图使用，以Normal为单位避免重复
## 按照Normal分类
lollipop_variant_sample <- unique(
    lollipop_variant[, c( "ID" , "gene" , "refseq" , "chromosome" , "start" , "REF" , "ALT" , "aachange" , "class" , "Class" , "From" , "Molecular.subtype" ) ]
    )

## Tumor
out_name <- paste0(out_path , "/" , gene , ".UniqueNormal.IGC_DGC.tsv" )
write.table(lollipop_variant_sample , out_name , row.names=F , sep ="\t" , quote = F)

## IGC
out_name <- paste0(out_path , "/" , gene , ".IGC.tsv" )
write.table( subset(lollipop_variant_sample , Class == "IGC")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

## DGC
out_name <- paste0(out_path , "/" , gene , ".DGC.tsv" )
write.table( subset(lollipop_variant_sample , Class == "DGC")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

## GC
out_name <- paste0(out_path , "/" , gene , ".GC.tsv" )
write.table( subset(lollipop_variant_sample , Class %in% c("IGC" , "DGC"))[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

############################################################################################################
## CIN和GS的IGC和DGC

out_name <- paste0(out_path , "/" , gene , ".IGC.CIN.tsv" )
write.table( subset(lollipop_variant_sample , Class == "IGC" & Molecular.subtype == "CIN")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

out_name <- paste0(out_path , "/" , gene , ".IGC.GS.tsv" )
write.table( subset(lollipop_variant_sample , Class == "IGC" & Molecular.subtype == "GS")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

out_name <- paste0(out_path , "/" , gene , ".DGC.CIN.tsv" )
write.table( subset(lollipop_variant_sample , Class == "DGC" & Molecular.subtype == "CIN")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)

out_name <- paste0(out_path , "/" , gene , ".DGC.GS.tsv" )
write.table( subset(lollipop_variant_sample , Class == "DGC" & Molecular.subtype == "GS")[,c(1:9,11:12)] , out_name , row.names=F , sep ="\t" , quote = F)
